package kr.ac.hanyang.cs.gradPrj;

public class LinearPredictiveCoding {

    double lambda = 0.0;

    int numDimensions = 10;

    /**
     * Basic constructor for LPC that sets definition, dependencies, and offsets
     * field.
     */
    public LinearPredictiveCoding() {
    }

    /**
     * Code taken from www.musicdsp.org.
     * <p>
     * mail.mutagene.net.2005. <i>Music dsp source archive</i> [online] [cited
     * May 10, 2005] Available on world wide web
     * (http://musicdsp.org/archive.php?classid=2#137)
     * 
     * @param samples
     *            The samples to extract the feature from.
     * @param sampling_rate
     *            The sampling rate that the samples are encoded with.
     * @return The extracted feature value(s).
     * @throws Exception
     *             Throws an informative exception if the feature cannot be
     *             calculated.
     * @see jAudioFeatureExtractor.AudioFeatures.FeatureExtractor#extractFeature(double[],
     *      double, double[][])
     */
    public double[] findCoefficient(short[] samples, double sampling_rate) {
            // find the order-P autocorrelation array, R, for the sequence x of
            // length L and warping of lambda
            // wAutocorrelate(&pfSrc[stIndex],siglen,R,P,0);

            double[] R = new double[numDimensions + 1];
            double K[] = new double[numDimensions];
            double A[] = new double[numDimensions];
            double[] dl = new double[samples.length];
            double[] Rt = new double[samples.length];
            double r1, r2, r1t;
            R[0] = 0;
            Rt[0] = 0;
            r1 = 0;
            r2 = 0;
            r1t = 0;
            for (int k = 0; k < samples.length; k++) {
                    Rt[0] += samples[k] * samples[k];

                    dl[k] = r1 - lambda * (samples[k] - r2);
                    r1 = samples[k];
                    r2 = dl[k];
            }
            for (int i = 1; i < R.length; i++) {
                    Rt[i] = 0;
                    r1 = 0;
                    r2 = 0;
                    for (int k = 0; k < samples.length; k++) {
                            Rt[i] += dl[k] * samples[k];

                            r1t = dl[k];
                            dl[k] = r1 - lambda * (r1t - r2);
                            r1 = r1t;
                            r2 = dl[k];
                    }
            }
            for (int i = 0; i < R.length; i++)
                    R[i] = Rt[i];

            // LevinsonRecursion(unsigned int P, float *R, float *A, float *K)
            double Am1[] = new double[62];

            if (R[0] == 0.0) {
                    for (int i = 1; i < numDimensions; i++) {
                            K[i] = 0.0;
                            A[i] = 0.0;
                    }
            } else {
                    double km, Em1, Em;
                    int k, s, m;
                    for (k = 0; k < numDimensions; k++) {
                            A[0] = 0;
                            Am1[0] = 0;
                    }
                    A[0] = 1;
                    Am1[0] = 1;
                    km = 0;
                    Em1 = R[0];
                    for (m = 1; m < numDimensions; m++) // m=2:N+1
                    {
                            double err = 0.0f; // err = 0;
                            for (k = 1; k <= m - 1; k++)
                                    // for k=2:m-1
                                    err += Am1[k] * R[m - k]; // err = err + am1(k)*R(m-k+1);
                            km = (R[m] - err) / Em1; // km=(R(m)-err)/Em1;
                            K[m - 1] = -km;
                            A[m] = km; // am(m)=km;
                            for (k = 1; k <= m - 1; k++)
                                    // for k=2:m-1
                                    A[k] = Am1[k] - km * Am1[m - k]; // am(k)=am1(k)-km*am1(m-k+1);
                            Em = (1 - km * km) * Em1; // Em=(1-km*km)*Em1;
                            for (s = 0; s < numDimensions; s++)
                                    // for s=1:N+1
                                    Am1[s] = A[s]; // am1(s) = am(s)
                            Em1 = Em; // Em1 = Em;
                    }
            }
            return K;
    }

}
